l05-06

library(leaflet)
popup = c("Robin", "Jakub", "Jannes")
leaflet() %>%
  addProviderTiles("Stamen.Toner") %>%
  addCircles(lng = c(-3, 23, 11),
             lat = c(52, 53, 49), 
             popup = popup)

#chapter2

library(sf) 

## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

library(raster) 

## Loading required package: sp

library(spData) # load geographic data
library(spDataLarge) # load larger geographic data

world

## Simple feature collection with 177 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS:  WGS 84
## # A tibble: 177 x 11
##    iso_a2 name_long continent region_un subregion type  area_km2     pop lifeExp
##    <chr>  <chr>     <chr>     <chr>     <chr>     <chr>    <dbl>   <dbl>   <dbl>
##  1 FJ     Fiji      Oceania   Oceania   Melanesia Sove…   1.93e4  8.86e5    70.0
##  2 TZ     Tanzania  Africa    Africa    Eastern … Sove…   9.33e5  5.22e7    64.2
##  3 EH     Western … Africa    Africa    Northern… Inde…   9.63e4 NA         NA  
##  4 CA     Canada    North Am… Americas  Northern… Sove…   1.00e7  3.55e7    82.0
##  5 US     United S… North Am… Americas  Northern… Coun…   9.51e6  3.19e8    78.8
##  6 KZ     Kazakhst… Asia      Asia      Central … Sove…   2.73e6  1.73e7    71.6
##  7 UZ     Uzbekist… Asia      Asia      Central … Sove…   4.61e5  3.08e7    71.0
##  8 PG     Papua Ne… Oceania   Oceania   Melanesia Sove…   4.65e5  7.76e6    65.2
##  9 ID     Indonesia Asia      Asia      South-Ea… Sove…   1.82e6  2.55e8    68.9
## 10 AR     Argentina South Am… Americas  South Am… Sove…   2.78e6  4.30e7    76.3
## # … with 167 more rows, and 2 more variables: gdpPercap <dbl>,
## #   geom <MULTIPOLYGON [°]>

plot(world)

## Warning: plotting the first 9 out of 10 attributes; use max.plot = 10 to plot
## all

names(world)

##  [1] "iso_a2"    "name_long" "continent" "region_un" "subregion" "type"     
##  [7] "area_km2"  "pop"       "lifeExp"   "gdpPercap" "geom"

#world["lifeExp"]
world %>% dplyr::select(lifeExp) %>% st_drop_geometry()

## # A tibble: 177 x 1
##    lifeExp
##  *   <dbl>
##  1    70.0
##  2    64.2
##  3    NA  
##  4    82.0
##  5    78.8
##  6    71.6
##  7    71.0
##  8    65.2
##  9    68.9
## 10    76.3
## # … with 167 more rows

#st_drop_geometry(world["lifeExp"])

class(world)

## [1] "sf"         "tbl_df"     "tbl"        "data.frame"

plot(world["lifeExp"])

world_asia = world[world$continent == "Asia", ]
asia = st_union(world_asia)

## although coordinates are longitude/latitude, st_union assumes that they are planar

#> although coordinates are longitude/latitude, st_union assumes that they are planar
#world_asia
#asia

#plot(world_asia)
plot(asia)

plot(world["pop"],reset = FALSE)
plot(asia, add = TRUE, col="red")

multipoint_matrix = rbind(c(5, 2), c(1, 3), c(3, 4), c(3, 2))
multi_point <- st_multipoint(multipoint_matrix)

multipoint_matrix

##      [,1] [,2]
## [1,]    5    2
## [2,]    1    3
## [3,]    3    4
## [4,]    3    2

class(multi_point)

## [1] "XY"         "MULTIPOINT" "sfg"

lnd_point = st_point(c(0.1, 51.5))                 # sfg object
class(lnd_point)

## [1] "XY"    "POINT" "sfg"

lnd_geom = st_sfc(lnd_point, crs = 4326)           # sfc object
class(lnd_geom)

## [1] "sfc_POINT" "sfc"

lnd_attrib = data.frame(                           # data.frame object
  name = "London",
  temperature = 25,
  date = as.Date("2017-06-21")
  )
lnd_attrib

##     name temperature       date
## 1 London          25 2017-06-21

lnd_sf = st_sf(lnd_attrib, geometry = lnd_geom)    # sf object
lnd_sf

## Simple feature collection with 1 feature and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 0.1 ymin: 51.5 xmax: 0.1 ymax: 51.5
## Geodetic CRS:  WGS 84
##     name temperature       date         geometry
## 1 London          25 2017-06-21 POINT (0.1 51.5)

#load data
raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge")
new_raster = raster(raster_filepath)
new_raster

## class      : RasterLayer 
## dimensions : 457, 465, 212505  (nrow, ncol, ncell)
## resolution : 0.0008333333, 0.0008333333  (x, y)
## extent     : -113.2396, -112.8521, 37.13208, 37.51292  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : /usr/local/lib/R/site-library/spDataLarge/raster/srtm.tif 
## names      : srtm 
## values     : 1024, 2892  (min, max)

plot(new_raster)

multi_raster_file = system.file("raster/landsat.tif", package = "spDataLarge")
multi_raster_file

## [1] "/usr/local/lib/R/site-library/spDataLarge/raster/landsat.tif"

r_brick = brick(multi_raster_file)
r_brick

## class      : RasterBrick 
## dimensions : 1428, 1128, 1610784, 4  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 301905, 335745, 4111245, 4154085  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs 
## source     : /usr/local/lib/R/site-library/spDataLarge/raster/landsat.tif 
## names      : landsat.1, landsat.2, landsat.3, landsat.4 
## min values :      7550,      6404,      5678,      5252 
## max values :     19071,     22051,     25780,     31961

plot(r_brick)

plotRGB(r_brick[[1:3]]/max(r_brick[[1:3]]),3,2,1,scale=T)

#This is a dataset containing the four bands (2, 3, 4, 5) of the Landsat 8 image for the area of Zion National Park

ndvi = (r_brick[[4]] - r_brick[[3]])/(r_brick[[3]]+r_brick[[4]])
plot(ndvi)

#chapter3

library(sf)
library(raster)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(stringr) # for working with strings (pattern matching)
library(tidyr)   # for unite() and separate()

## 
## Attaching package: 'tidyr'

## The following object is masked from 'package:raster':
## 
##     extract

sel_area = world$area_km2 < 10000
small_countries = world[sel_area, ]
small_countries

## Simple feature collection with 7 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -67.24243 ymin: -16.59785 xmax: 167.8449 ymax: 50.12805
## Geodetic CRS:  WGS 84
## # A tibble: 7 x 11
##   iso_a2 name_long  continent region_un subregion type  area_km2     pop lifeExp
##   <chr>  <chr>      <chr>     <chr>     <chr>     <chr>    <dbl>   <dbl>   <dbl>
## 1 PR     Puerto Ri… North Am… Americas  Caribbean Depe…    9225. 3534874    79.4
## 2 PS     Palestine  Asia      Asia      Western … Disp…    5037. 4294682    73.1
## 3 VU     Vanuatu    Oceania   Oceania   Melanesia Sove…    7490.  258850    71.7
## 4 LU     Luxembourg Europe    Europe    Western … Sove…    2417.  556319    82.2
## 5 <NA>   Northern … Asia      Asia      Western … Sove…    3786.      NA    NA  
## 6 CY     Cyprus     Asia      Asia      Western … Sove…    6207. 1152309    80.2
## 7 TT     Trinidad … North Am… Americas  Caribbean Sove…    7738. 1354493    70.4
## # … with 2 more variables: gdpPercap <dbl>, geom <MULTIPOLYGON [°]>

plot(small_countries["pop"])

world %>% filter(area_km2 < 10000) %>% 
  dplyr::select(pop) %>% 
  plot()

##1.GDPperCapitaがtop10の国をピックアップして地図化してみましょう
#top_n
world %>% top_n(n=10,wt=gdpPercap) %>% 
  dplyr::select(gdpPercap) %>% 
  plot()

##2.GDPperCapitaがworst10の国をピックアップして地図化してみましょう
#top_n
world %>% top_n(n=-10,wt=gdpPercap) %>% 
  dplyr::select(gdpPercap) %>% 
  plot()

##3.日本はGDPperCapitaで上から何番目でしょう
world %>% dplyr::select(name_long,gdpPercap)

## Simple feature collection with 177 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS:  WGS 84
## # A tibble: 177 x 3
##    name_long     gdpPercap                                                  geom
##    <chr>             <dbl>                                    <MULTIPOLYGON [°]>
##  1 Fiji              8222. (((180 -16.06713, 180 -16.55522, 179.3641 -16.80135,…
##  2 Tanzania          2402. (((33.90371 -0.95, 34.07262 -1.05982, 37.69869 -3.09…
##  3 Western Saha…       NA  (((-8.66559 27.65643, -8.665124 27.58948, -8.6844 27…
##  4 Canada           43079. (((-122.84 49, -122.9742 49.00254, -124.9102 49.9845…
##  5 United States    51922. (((-122.84 49, -120 49, -117.0312 49, -116.0482 49, …
##  6 Kazakhstan       23587. (((87.35997 49.21498, 86.59878 48.54918, 85.76823 48…
##  7 Uzbekistan        5371. (((55.96819 41.30864, 55.92892 44.99586, 58.50313 45…
##  8 Papua New Gu…     3709. (((141.0002 -2.600151, 142.7352 -3.289153, 144.584 -…
##  9 Indonesia        10003. (((141.0002 -2.600151, 141.0171 -5.859022, 141.0339 …
## 10 Argentina        18798. (((-68.63401 -52.63637, -68.25 -53.1, -67.75 -53.85,…
## # … with 167 more rows

#gdpPercapをソートしたい
#sort(world$gdpPercap,decreasing = T)
gdp_sort <- world %>% arrange(desc(gdpPercap)) %>% 
  dplyr::select(gdpPercap) %>% 
  st_drop_geometry()

#日本のgdpPercapを取り出したい
jp_gdp <- world %>% filter(name_long == "Japan") %>% 
  dplyr::select(gdpPercap) %>% 
  st_drop_geometry() %>% 
  as.numeric()

#gdpPercapが上から何番目か調べられれば良い
#which()
which(gdp_sort==jp_gdp)

## [1] 22